% title:             insurance model with decreased variation of costs
% author:            neale mahoney
% lastest revision:  7/03/14


%% Setup

clear all; clc;

% cd(<folder directory>);

version = 'v4';
%% Primatives

fprintf('Displaying primatives');

n_i = 16e3;             %number of individuals
n_c = 1e3;              %number of cost shocks per individual
smooth_type = 'lowess'; %moving, lowess, rlowess
bw = .1;                %bandsmooth for moving average smoother

%% Joint distribution of risk aversion and expected medical costs

fprintf('\n\nBaseline variation of costs\n')

% %distribution of absoluate risk aversion from Table 3 of Handel, Hendel and Whinston (2013)
m = 4.39e-4;
v = (6.63e-5)^2;
mu_alpha = log((m^2)/sqrt(v+m^2));
sigma_alpha = sqrt(log(v/(m^2)+1));

%gamble interpretation
f = zeros(100,1);
for i = 1:100,
    f(i) = cara_function(0,m,[-100 i],[0 0]);
end

[value index] = min(abs(f));
fprintf('Mean risk adversion makes individual different between 50-50 gamble for (100,%i) and 0 \n', index)
  
% distribution of health type
m = .9794951;
v = 1.377593^2;
mu_l = log((m^2)/sqrt(v+m^2));
sigma_l = sqrt(log(v/(m^2)+1));

%joint normal
mu = [mu_alpha; mu_l];
rho = 0;
cov = rho*sigma_alpha*sigma_l;
sigma = [sigma_alpha^2 cov; cov sigma_l^2];
X = mvnrnd(mu,sigma,n_i);

%distribution of risk preferences and health type
alpha = exp(X(:,1));
lambda = exp(X(:,2));

%distribution of realized costs
m = 3138.765; 
v = 10125.74^2;
mu_c = log((m^2)/sqrt(v+m^2));
sigma_c = sqrt(log(v/(m^2)+1));
rho = .498481; 

%distribution of realized cost conditional on health type
c = zeros(n_i,n_c);
part_1 = mu_c;
part_2 = rho*sigma_c*(X(:,2)-mu_l)/sigma_l;

for i = 1:n_i   
    part_3 = sqrt((1-rho^2))*sigma_c*randn(1,n_c);
    c(i,:) = exp(part_1 + part_2(i) + part_3);
end

%av 60 plan: low quality plan that covers 60% of medical cost on average
c_oop_60 = min(min(c,500+.4*(c-500)),8e3);  %500 deductible, 8000 out of pocket maximum
c_ins_60 = c - c_oop_60;
c_ins_mean_60 = mean(c_ins_60,2);

fprintf('AV 60 mean plan cost is %i \n',round(mean(c_ins_mean_60)))
fprintf('AV 60 mean oop cost is %i \n',round(mean(reshape(c_oop_60,n_i*n_c,1))))
fprintf('AV 60 AV is %f \n',mean(c_ins_mean_60)/mean(reshape(c,n_i*n_c,1)))

% av 90 plan: high quality plan that covers 90% of medical cost on average
c_oop_90 = min(.1*c,8e3);   %no deductible, 8000 out of pocket maximum
c_ins_90 = c - c_oop_90;
c_ins_mean_90 = mean(c_ins_90,2);

fprintf('AV 90 mean plan cost is %i \n',round(mean(c_ins_mean_90)))
fprintf('AV 90 mean oop cost is %i \n',round(mean(reshape(c_oop_90,n_i*n_c,1))))
fprintf('AV 90 AV is %f \n',mean(c_ins_mean_90)/mean(reshape(c,n_i*n_c,1)))
fprintf('Delta AC is %f \n', round(mean(c_ins_mean_90)-mean(c_ins_mean_60)))

%allow bankruptcy to cap out of pocket risk
c = min(c,20e3);


%% Decrease variation of costs
fprintf('\n\nDecreased variation of costs\n')

m = 2210;
v = (10125.74*.2)^2;
mu_c_scale = log((m^2)/sqrt(v+m^2));
sigma_c_scale = sqrt(log(v/(m^2)+1));
rho = .498481;

c_scale = zeros(n_i,n_c);
part_1 = mu_c_scale;
part_2 = rho*sigma_c_scale*(X(:,2)-mu_l)/sigma_l;

for i = 1:n_i
    part_3 = sqrt((1-rho^2))*sigma_c_scale*randn(1,n_c);
    c_scale(i,:) = exp(part_1 + part_2(i) + part_3);
end

%av 60 plan
c_oop_60_scale = min(min(c_scale,500+.4*(c_scale-500)),8e3);
c_ins_60_scale = c_scale - c_oop_60_scale;
c_ins_mean_60_scale = mean(c_ins_60_scale,2);

fprintf('AV 60 mean plan cost is %i \n',round(mean(c_ins_mean_60_scale)))
fprintf('AV 60 mean oop cost is %i \n',round(mean(reshape(c_oop_60_scale,n_i*n_c,1))))
fprintf('AV 60 AV is %f \n',mean(c_ins_mean_60)/mean(reshape(c_scale,n_i*n_c,1)))

% av 90 plan
c_oop_90_scale = min(.1*c_scale,8e3);
c_ins_90_scale = c_scale - c_oop_90_scale;
c_ins_mean_90_scale = mean(c_ins_90_scale,2);

fprintf('AV 90 mean plan cost is %i \n',round(mean(c_ins_mean_90_scale)))
fprintf('AV 90 mean oop cost is %i \n',round(mean(reshape(c_oop_90_scale,n_i*n_c,1))))
fprintf('AV 90 AV is %f \n',mean(c_ins_mean_90_scale)/mean(reshape(c_scale,n_i*n_c,1)))
fprintf('Delta AC is %f \n', round(mean(c_ins_mean_90_scale)-mean(c_ins_mean_60_scale)))

c_scale = min(c_scale,20e3);

%% wtp
fprintf('\n\nEstimating WTP\n\n');

v_grid = 0:10:5e4;
v = zeros(n_i,1);
for i = 1:n_i,
    %v: willingness to pay additional premium such that 
    %the customer is indifferent between actual costs of high and low quality plan
    if mod(i,25) == 0, fprintf('Estimating WTP for observation %i \n',i); end
    v(i) = wtp_function(alpha(i),v_grid,c_oop_90(i,:),c_oop_60(i,:));

end

%% Equilibrium - Baseline

%sort matrix X in descending order by wtp and smooth
wtp_shift = 900;
X = [v + wtp_shift c_ins_mean_90-c_ins_mean_60];
X = sortrows(X,-1);

v_sort = smooth(X(:,1),bw,smooth_type);
mc = smooth(X(:,2),bw,smooth_type);
q = transpose((1:n_i)./n_i);

plot_switch = 1;
plot_number = 1;
plot_dir = '';
plot_name = 'correlations_baseline';

theta = .5;
out_baseline = equilibrium_function_v2(...
    v_sort,mc,q,theta,...
    bw,smooth_type,plot_switch,plot_number,plot_dir,plot_name,version); 

%% Equilibrium - Decreased variation of costs
%sort matrix X in descending order by wtp and smooth
wtp_shift = 900;
X = [v + wtp_shift c_ins_mean_90_scale-c_ins_mean_60_scale];
X = sortrows(X,-1);

v_sort = smooth(X(:,1),bw,smooth_type);
mc_scale = smooth(X(:,2),bw,smooth_type);
q = transpose((1:n_i)./n_i);

plot_switch = 1;
plot_number = 2;
plot_dir = '';
plot_name = 'correlations_alt';

%market power: theta = 0.5
%this corresponds to Cournot competition with two high-quality plans
theta = .5;
out_alt = equilibrium_function_v2(...
    v_sort,mc_scale,q,theta,...
    bw,smooth_type,plot_switch,plot_number,plot_dir,plot_name,version);







